library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
#library(corrplot)
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
source("~/Documents/GitHub/RISKPLOTS/CODE/auxplots.R")
data(cgd)
data <- model.frame(Surv(time,status)~.*.,lung,na.action=NULL)
colnames(data) <-str_replace_all(colnames(data),":","_")
colnames(data) <-str_replace_all(colnames(data),"\\.","_")
data$inst <- NULL
data$`Surv(time, status)` <- NULL
dataLung <- cbind(time=lung$time/365,status=lung$status-1,data)
dataLung <- dataLung[complete.cases(dataLung),]
pander::pander(table(dataLung$status))
pander::pander(summary(dataLung$time))
| 0.0137 |
0.4788 |
0.7356 |
0.8494 |
1.14 |
2.8 |
Modeling
ml <- BSWiMS.model(Surv(time,status)~1,data=dataLung,NumberofRepeats = 10)
[++++++++++++++++++++++++++]..
sm <- summary(ml)
pander::pander(sm$coefficients)
Table continues below
| ph_ecog |
0.4323 |
1.194 |
1.541 |
1.988 |
0.6786 |
| sex |
-0.4591 |
0.456 |
0.6319 |
0.8756 |
0.6488 |
| pat_karno |
-0.001767 |
0.9968 |
0.9982 |
0.9997 |
0.506 |
| ph_karno |
-3.482e-07 |
1 |
1 |
1 |
0.5774 |
Table continues below
| ph_ecog |
0.6488 |
0.6012 |
0.6012 |
0.6196 |
0.5995 |
| sex |
0.6786 |
0.6012 |
0.6196 |
0.6012 |
0.5995 |
| pat_karno |
0.7202 |
0.506 |
0.5855 |
0.5 |
0.5855 |
| ph_karno |
0.7202 |
0.5774 |
0.57 |
0.5 |
0.57 |
| ph_ecog |
0.04491 |
0.4048 |
3.326 |
2.485 |
-0.02005 |
1 |
| sex |
0.0285 |
0.4783 |
2.758 |
2.85 |
-0.00167 |
1 |
| pat_karno |
0.02917 |
0.3418 |
2.439 |
2.243 |
0.08546 |
1 |
| ph_karno |
0.01425 |
0.2799 |
2.221 |
1.642 |
0.06998 |
0.6 |
Cox Model Performance
Here we evaluate the model using the RRPlot() function.
The evaluation of the raw Cox model with RRPlot()
Here we will use the predicted event probability assuming a baseline
hazard for events withing 5 years
index <- predict(ml,dataLung)
timeinterval <- 1 # One year
h0 <- sum(dataLung$status & dataLung$time < timeinterval)
h0 <- h0/nrow(subset(dataLung,time > timeinterval | status==1))
rdata <- cbind(dataLung$status,ppoisGzero(index,h0))
rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=dataLung$time,
title="Raw Train: Lung Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
As we can see the Observed probability as well as the Time vs. Events
are not calibrated.
Cox Calibration
op <- par(no.readonly = TRUE)
calprob <- CoxRiskCalibration(ml,dataLung,"status","time")
( 2.32615 , 781.9918 , 1.468558 , 121 , 140.6966 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
Crossvalidation
rcv <- randomCV(theData=dataLung,
theOutcome = Surv(time,status)~1,
fittingFunction=BSWiMS.model,
trainFraction = 0.9,
repetitions=200,
classSamplingType = "LOO"
)
## .[+++].[++].[++].[++].[+++].[++-].[++].[+++].[+++].[++++]10 Tested: 20 Avg. Selected: 3.6 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.4271312
##
.[++++].[++++].[+++].[++].[++++].[+++].[+++].[+].[++].[++]20 Tested: 40 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.422639
##
.[+++].[+++].[++].[++].[+++].[+++].[+++].[+++].[+++].[++]30 Tested: 60 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.4261273
##
.[++-].[++].[++].[++].[++].[+++].[+++].[++].[++].[+++]40 Tested: 80 Avg. Selected: 3.6 Min Tests: 1 Max Tests: 1 Mean Tests: 1 . MAD: 0.4419494
##
.[+++].[+++].[+++].[+++].[++-].[+++].[+++].[+++].[+++].[+++]50 Tested: 97 Avg. Selected: 3.66 Min Tests: 1 Max Tests: 2 Mean Tests: 1.030928 . MAD: 0.4486325
##
.[+++].[+++].[++-].[++].[+++].[+++].[++++].[+++].[+++].[+++]60 Tested: 107 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 2 Mean Tests: 1.121495 . MAD: 0.4708602
##
.[++].[+++].[+++].[++].[++].[+++].[++].[+++].[+++].[++-]70 Tested: 117 Avg. Selected: 3.671429 Min Tests: 1 Max Tests: 2 Mean Tests: 1.196581 . MAD: 0.4750843
##
.[+].[+++].[+++].[+++].[+++].[++].[+++].[++].[+++].[++]80 Tested: 127 Avg. Selected: 3.65 Min Tests: 1 Max Tests: 2 Mean Tests: 1.259843 . MAD: 0.4704381
##
.[++].[+++].[++].[+++].[+++].[+++].[+++].[+++].[+++].[+++]90 Tested: 137 Avg. Selected: 3.666667 Min Tests: 1 Max Tests: 2 Mean Tests: 1.313869 . MAD: 0.4791172
##
.[+++].[+++].[+++].[+++].[+++].[+++].[++].[+++].[+++].[++]100 Tested: 147 Avg. Selected: 3.68 Min Tests: 1 Max Tests: 3 Mean Tests: 1.360544 . MAD: 0.4774504
##
.[++].[+++].[+++].[++++].[++++].[+++].[+++].[++].[++++].[+++]110 Tested: 157 Avg. Selected: 3.718182 Min Tests: 1 Max Tests: 3 Mean Tests: 1.401274 . MAD: 0.4859387
##
.[++].[+++].[++].[+++].[+++].[+++].[+++].[++].[++].[++]120 Tested: 167 Avg. Selected: 3.7 Min Tests: 1 Max Tests: 3 Mean Tests: 1.437126 . MAD: 0.4880085
##
.[+++].[++].[+++].[++].[++].[++-].[++++].[++].[++].[+++]130 Tested: 168 Avg. Selected: 3.676923 Min Tests: 1 Max Tests: 3 Mean Tests: 1.547619 . MAD: 0.4865085
##
.[+++].[+++].[+++].[+++].[+++].[+++].[+].[+++].[++-].[+++]140 Tested: 168 Avg. Selected: 3.678571 Min Tests: 1 Max Tests: 3 Mean Tests: 1.666667 . MAD: 0.4884544
##
.[+++].[+++].[+++].[++-].[+++].[+++].[++].[++].[+++].[++]150 Tested: 168 Avg. Selected: 3.673333 Min Tests: 1 Max Tests: 4 Mean Tests: 1.785714 . MAD: 0.4898956
##
.[++++].[+++].[+++].[+++].[++-].[+++].[+++].[++].[+++].[++-]160 Tested: 168 Avg. Selected: 3.68125 Min Tests: 1 Max Tests: 4 Mean Tests: 1.904762 . MAD: 0.4911541
##
.[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++].[++-].[++]170 Tested: 168 Avg. Selected: 3.688235 Min Tests: 1 Max Tests: 4 Mean Tests: 2.02381 . MAD: 0.4897295
##
.[+++].[+++].[+++].[++-].[+++].[+++].[+++].[++].[++].[+++]180 Tested: 168 Avg. Selected: 3.688889 Min Tests: 1 Max Tests: 4 Mean Tests: 2.142857 . MAD: 0.4894258
##
.[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++].[+++]190 Tested: 168 Avg. Selected: 3.705263 Min Tests: 1 Max Tests: 5 Mean Tests: 2.261905 . MAD: 0.4880488
##
.[++].[+++].[+++].[++-].[++].[+++].[+++].[+++].[++++].[+++]200 Tested: 168 Avg. Selected: 3.71 Min Tests: 1 Max Tests: 5 Mean Tests: 2.380952 . MAD: 0.4906368
##
bbx <- boxplot(rcv$survTestPredictions[,1]~rownames(rcv$survTestPredictions),plot=FALSE)
times <- bbx$stats[3,]
status <- boxplot(rcv$survTestPredictions[,2]~rownames(rcv$survTestPredictions),plot=FALSE)$stats[3,]
prob <- ppoisGzero(boxplot(rcv$survTestPredictions[,4]~rownames(rcv$survTestPredictions),plot=FALSE)$stats[3,],h0)
rdatacv <- cbind(status,prob)
rownames(rdatacv) <- bbx$names
rrAnalysisTest <- RRPlot(rdatacv,atRate=c(0.90),
timetoEvent=times,
title="Test: Lung Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
Calibrating the test results
rdatacv <- cbind(status,prob,times)
calprob <- CalibrationProbPoissonRisk(rdatacv)
## ( 2.356124 , 781.9918 , 1.487481 , 121 , 143.3503 )
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
timeinterval <- calprob$timeInterval;
rdata <- cbind(status,calprob$prob)
rrAnalysisTest <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=times,
title="Calibrated Test: Lung",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
### Calibrated Test Performance
pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
Threshold values
| Thr |
0.8618 |
0.6585 |
0.5694 |
0.5147 |
0.5147 |
| RR |
1.186 |
2.958 |
3.293 |
14.58 |
14.58 |
| RR_LCI |
0.972 |
1.387 |
1.392 |
0.03466 |
0.03466 |
| RR_UCI |
1.447 |
6.309 |
7.792 |
6133 |
6133 |
| SEN |
0.1983 |
0.9587 |
0.9669 |
1 |
1 |
| SPE |
0.8936 |
0.2979 |
0.2766 |
0.04255 |
0.04255 |
| BACC |
0.546 |
0.6283 |
0.6218 |
0.5213 |
0.5213 |
| NetBenefit |
-0.04266 |
0.3119 |
0.4289 |
0.4362 |
0.4362 |
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
| 1.251 |
1.038 |
1.494 |
0.01672 |
pander::pander(rrAnalysisTest$OERatio$atThrEstimates,caption="O/E Ratio")
O/E Ratio
| Total |
121 |
100.4 |
144.6 |
96.74 |
1.251 |
1.038 |
1.494 |
0.01672 |
| low |
98 |
79.56 |
119.4 |
73.01 |
1.342 |
1.09 |
1.636 |
0.00489 |
| 90% |
23 |
14.58 |
34.51 |
24.12 |
0.9537 |
0.6046 |
1.431 |
0.9189 |
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Mean")
O/E Mean
| 0.9892 |
0.9892 |
0.9567 |
1.022 |
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
| 0.9732 |
0.9733 |
0.9627 |
0.9828 |
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
| 0.6027 |
0.6017 |
0.5349 |
0.6684 |
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
| 0.6052 |
0.5027 |
0.7077 |
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
| 0.1901 |
0.1245 |
0.2714 |
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
| 0.8936 |
0.769 |
0.9645 |
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
| 0.8632 |
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 2.346444 on 1 degrees of freedom, p =
0.125569
| class=0 |
140 |
98 |
103.8 |
0.3276 |
2.346 |
| class=1 |
28 |
23 |
17.17 |
1.981 |
2.346 |
